ESM 244 Lab 1

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library(tidyverse)
library(here)
library(sf)
library(tmap)

Read in the data

sf_trees <- read_csv(here('data', 'sf_trees', 'sf_trees.csv'))
# read.csv is the built in function; read_csv is an improved version from the tidyverse, has extra functionality and is faster
# summary(sf_trees); names(sf_trees)
# sf_trees %>% group_by(legal_status) %>% summarize(tree_count=n())
top_5_status <- sf_trees %>% 
  group_by(legal_status) %>% 
  summarize(tree_count = n()) %>% 
  slice_max(tree_count, n = 5) %>% 
  arrange(-tree_count)

ggplot(data = top_5_status, aes(x = fct_reorder(legal_status, tree_count),
                                y = tree_count)) +
  geom_col() +
  labs(y = 'Tree count', x = 'Legal Status') +
  coord_flip() +
  theme_minimal()

# ggplot is more general than ggpubr 

Keep observations where legal status is ‘Permitted Site’ and caretaker is ‘MTA’

permitted_mta <- sf_trees %>% 
  filter(legal_status == "Permitted Site" & caretaker == "MTA")
# == means matches. single equals sign assigns a value

Keep Blackwood acacia trees, only legal_status, latitute, longitude columns

blackwood_acacia <- sf_trees %>% 
  filter(str_detect(species, "Blackwood Acacia")) %>% 
  select(legal_status, date, latitude, longitude)

### make a plot of SF acacia
ggplot(data = blackwood_acacia, aes(x = longitude, y = latitude)) +
  geom_point()

### Use tidyr::separate() to separate one column into multiple, and tidyr::unite() to rejoin

sf_trees_sep <- sf_trees %>% 
  separate(species, into = c('spp_scientific', 'spp_common'), sep = ' :: ')

sf_trees_unite <- sf_trees %>% 
  unite("id_status", tree_id:species, sep = '_NEW_')

Make some actual maps !

blackwood_acacia_sf <- blackwood_acacia %>% 
  drop_na(longitude, latitude) %>% 
  st_as_sf(coords = c('longitude', 'latitude'))
# sf stands for spatial features
st_crs(blackwood_acacia_sf) <- 4326 # coordinate reference system indicating basic latitude and longitude
ggplot(data=blackwood_acacia_sf) +
  geom_sf(color = 'darkgreen') + # always looking for a geometry column and will use it as aes
  theme_minimal()

read in SF streets data

sf_map_sf <- read_sf(here('data', 'sf_map', 'tl_2017_06075_roads.shp')) %>% 
  st_transform(4326) # change the coordinate reference system
st_crs(sf_map_sf) # check to make sure it has changed

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

ggplot() +
  geom_sf(data = sf_map_sf, size = 0.1, color = 'darkgrey') +
  geom_sf(data = blackwood_acacia_sf, color = 'red', size = 0.5) +
  theme_void()

  labs(title = 'Blackwood acacias in San Francisco')

## $title
## [1] "Blackwood acacias in San Francisco"
## 
## attr(,"class")
## [1] "labels"

Interactive map!

tmap_mode('view')
tm_shape(blackwood_acacia_sf) +
  tm_dots()